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Abstract 

We give a brief account of the development of methods to include thermal 
fluctuations into lattice Boltzmann algorithms. Emphasis is put on our recent 
work (Phys. Rev. E 76, 036704 (2007)) which provides a clear understanding 
in terms of statistical mechanics. 
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The lattice Boltzmann (LB) equation has, in the last few decades, 
emerged as a powerful tool to solve fluid dynamics problems numerically 
[ll, [ij. The algorithm is a fully discretized version of the Boltzmann equa- 
tion, known from the kinetic theory of gases. Space r is discretized in terms 
of a regular (usually simple-cubic) lattice with spacing 6, time t in terms of a 
time step h, and velocity space in terms of a small set of velocities q that are 
chosen such that q/i is a vector which connects two nearby lattice sites. For 
example, the popular D3Q19 model [sl employs nineteen velocities, corre- 
sponding to zero and the six nearest and twelve next-nearest neighbors on a 
simple-cubic lattice. The central quantities on which the algorithm operates 
are the populations nj(r, t), representing the mass density corresponding to 
velocity q, such that the total mass density p{r,t) at the site r at time t is 
given by 

p(f,t) = ^n,(f,t). (1) 
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Similarly, the momentum density is obtained as the first velocity moment, 



j{f,t) = Y^n,{f,t)ci, (2) 

i 

and the hydrodynamic flow velocity is given by 

u{r,t) = ——. 3 

The algorithm is then described by the lattice Boltzmann equation 

ni{f+ Cih, t + h)= n*{f, t) = ni{f, t) + {{ni{f, t)}) . (4) 

The collision operator Aj modifies the populations on the site {{rii} denotes 
the set of all populations on the site), such that mass and momentum are 
conserved. Energy conservation is not taken into account, since we are here 
interested in an isothermal version, where the temperature instead of the 
energy is fixed (formally, this corresponds to a system with infinite heat 
conductivity). The conservation equations therefore read 

J2A, = Y,^^C^ = 0. (5) 

i i 

This results in a set of post-coUisional populations n*, which are then prop- 
agated to the neighboring sites. 

In most applications, it is assumed that Aj is a deterministic variable, 
i. e. that it can be calculated in a unique fashion from the populations 
Hi (r , t) . This is very much in spirit of the original continuum Boltzmann 
equation, and applicable to many practical problems of fiuid flow. However, 
for soft-matter applications, where one is interested in Brownian motion of 
suspended particles, or similar phenomena, this is not sufficient. Rather, one 
must take into account that here both the lattice spacing b and the time 
step h are so small that on these scales thermal fluctuations are sizeable and 
cannot be viewed as just averaged out. Indeed, assuming that the underlying 
physical model is an ideal gas, one can see this rather easily by starting from 
the equation of state 

ksT = mpcl, (6) 

where is Boltzmann's constant, T is the absolute temperature, rrip is 
the mass of a gas particle, and Cg is the isothermal speed of sound (c^ — 
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p/ Pi where p is the thermodynamic pressure). Usually, is chosen as an 
adjustable parameter, picked in such a way that — even in nonequilibrium 
situations like shear flow — the typical flow velocity u is small compared 
to Cs- This is the condition of low Mach number flow, which is needed 
because of the restricted velocity space (note that Cg is of the order of the 
Cj). Furthermore, the physics of the problem usually dictates the values of 
ksT and p — for example, we may assume that we study water at room 
temperature. Equation [6] then allows us to determine the mass of a gas 
particle, m^, which, in turn, determines the number of particles on a lattice 
site (assuming a simple-cubic lattice in three dimensions), 

iV, = (7) 

If this number is very large, fluctuations will strongly average out, i. e. one 
can consider the single lattice site as a thermodynamic system. This is the 
case for typical engineering applications. However, if Np is comparable to 
unity, as it is the case for many soft-matter applications, then fluctuations 
are important, and must be taken into account in the algorithm. Since the 
system is an ideal gas, Np is a random variable whose probability distribution 
is Poisson. For such a distribution, the variance is identical to the mean, i. e. 
the relative importance of fluctuations is given by 



(we coined the word "Boltzmann number" for this parameter). We thus see 
that the degree of fluctuations is controlled by the degree of coarse-graining, 
through the lattice spacing b. It is also useful to introduce the parameter 

which may be called the thermal mass density. 

The question of how to actually implement these fluctuations in the col- 
lision operator Aj has found different answers during the last flfteen years, 
with increasing level of reflnement and understanding. In what follows, we 
wish to briefly outline these developments. Since all the material has been 
published previously, we would like to be brief, and refer the interested reader 
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to the original papers j3, [s], 0, 0] as well as to a recent review [sl , in which 
all the technical details have been worked out and explained in depth. 

The first implementation of a fluctuating lattice Boltzmann equation was 
by Ladd 0, He started from the well- understood deterministic version 
{Bo = 0), and added a stochastic term A- to the collision operator, with the 
requirement that this is consistent, on the macroscopic scale, with fluctuating 
hydrodynamics, as given by Landau and Lifshitz [9,]. 

Let us first discuss the deterministic version in some more detail. It is 
based upon a linearized collision operator, 

A, = ^L,,(n,-nf), (10) 
j 

where the matrix Lij contains constant elements, and is implicitly given via 
a diagonal representation (see below), while n^'' is the lattice analog to a 
velocity-dependent Maxwell-Boltzmann distribution: 

ip, u) = a-p ( 1 + ^ + - ^ 1 . (11) 



c 



Here Cs is the speed of sound, and the weights a'^^ > are normalized such 
that a'^' = 1. This notation has been chosen in order to emphasize that, 
for symmetry reasons, the weights only depend on the absolute values of the 
speeds q, but not on their direction. Furthermore, the weights are adjusted 
in such a way that n^'' satisfies the properties 



«r = p. (12) 

n'i'Si = j, (13) 

nfci^Ci = pc1l+pu0u=Il . (14) 



For D3Q19, this implies a'^^ = 1/3 for the rest population, a^* = 1/18 for the 
nearest neighbors, and a^* = 1/36 for the next-nearest neighbors. Further- 
more cl = (l/3)(6V/i2). 

Lij is implemented as follows: First, one transforms to so-called "modes", 
i. e. linear combinations of the rii which are adapted to the symmetry of the 
problem. The first ten modes have a direct hydrodynamic interpretation: 
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• Mode 0: Mass density p = Ylii'^i- 

• Modes 1-3: Momentum density ja = '^iUiCia; here a denotes a Carte- 
sian index. 

• Modes 4-9: Stresses Uais = 'Yli'^'^i^ia'^iP) which are conveniently decom- 
posed into trace and traceless part: = -|- ^^Sapli^^] here we use 
the Einstein summation convention. 

The additional modes (so-called "kinetic" or "ghost" modes) do not have a 
direct relation to hydrodynamics. In the D3Q19 model, there are nine such 
modes, which are explicitly listed in Ref. [8|. After having calculated the 
(pre-coUisional) modes, one leaves the conserved modes unchanged, while the 
other modes are linearly relaxed towards their local equilibrium value. The 
stresses are changed from pre- to post-coUisional values according to 

n:r = 7.n7, (15) 

where we use the notation n"^'' = Ui — rfj^ . The kinetic modes are defined in 
such a way that their equilibrium part is zero, and the action of Lij on them 
is, in the simplest version, just a projection, such that the post-collisional 
kinetic modes vanish. 

A Chapman- Enskog analysis shows that this procedure yields the Navier- 
Stokes equations of hydrodynamics in the limit of large length and time 
scales, with shear and bulk viscosities that are uniquely determined by the 
values of 7^ and 7;,, respectively. Linear stability requires |7s| < 1, |7fe| < 1, 
corresponding to positive values of the viscosities. 

This deterministic procedure was modified by Ladd jl,[5| by just changing 
Eq. da to 

n:r = 7.nr/ + ^a„ (le) 

'-'-aa ib'-'-aa ~ ^^aai 

with suitably chosen random stresses Rai3, while the treatment of the kinetic 
modes was left unchanged. The rationale behind this procedure was that 
the kinetic modes do not contribute to hydrodynamics, and the goal was 
to simulate the fluctuations correctly on the hydrodynamic scale. On this 
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scale, however, the fluctuating stresses Ra/3 that appear in the Navier-Stokes 
equation (different from Rap that appears in Eq. ^E^j satisfy the relations 



Raf^) = 0, (17) 

Ra(3 (r , t) R^s {f", t') ) = 2kBTr]a(3jS S (f - 6 {t - f) 

2kBT 

where rjafs-ys is the isotropic fourth-rank viscosity tensor, parameterized by 
shear and bulk viscosity, or the relaxation parameters 7^ and 7;,. In the last 
step, we have discretized the delta functions by the lattice parameter b and 
the time step h, as it is appropriate for a lattice simulation. 

One might expect that the LB noises are just given by Ra/3 = Rap- How- 
ever, this turns out not to be correct 0, 0]. Rather, the correct fluctuating 
LB stresses are obtained by a suitable modification of the amplitude. For 
the shear stresses one has, for example. 



{Rly) = {'^-ls?{Rly). (18) 

The same modification factor occurs for all other shear stresses, too, while 
the corresponding factor for the bulk stresses is (1 — 7f,)^. The reason has 
been explained in detail in Refs. jl, [sl; essentially the renormalization of 
the amplitude comes from the fact that Eq. [T7] describes the physics on 
a more coarse-grained time scale than Eq. [16] — the delta correlation in 
time is in LB replaced by an exponential decay. However, the time integral 
of the correlation functions must be the same in order to obtain the same 
macroscopic viscosities. 

Adhikari et al. 0] then generalized this procedure by not only thermal- 
izing the stresses, but also the kinetic modes, which were treated in a rather 
similar fashion to Eq. [161 The argument was that the relaxation of ki- 
netic modes introduces an additional dissipative mechanism into the system, 
which should be balanced by a compensating Langevin noise. A projection 
should be viewed as the limit of such a relaxation, with relaxation parameter 
7 — > 0, such that the fluctuation-dissipation relation should hold in this case, 
too. While this argument makes intuitive sense, and led to a substantially 
improved representation of the fluctuations at short length scales j6| , the the- 
oretical foundation of this procedure remained somewhat obscure (at least 
to the present authors). 
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In a recent publication [3] we have been able to resolve these questions by 
developing a first-principles theory of thermal fluctuations in LB models. The 
starting point was the observation that for a discrete system the concept of a 
fluctuation-dissipation theorem should rather be replaced by the concept of 
detailed balance as it applies to Monte Carlo simulations |10]. In order to be 
able to check whether an update rule satisfies or violates the detailed-balance 
condition, we therefore explicitly constructed the probability density for the 
random variables Ui on a site in thermal equilibrium. Taking advantage of 
the underlying picture of a gas of particles, we first transform from the Ui to 
variables Vi, the number of particles on the site which have velocity q (cf. 
Eqs. El and ED: 

y^ = -. (19) 
^^ 

In terms of these variables, the probability density (except for normalization, 
which is unimportant for our purposes) is written as 

p ({^^}) (n ^ (-'^*)) ^ ^"^^ -p^^i^ /^^^'^^ - ■ (20) 

The underlying picture is that of a "velocity bin" i in thermal contact with a 
huge reservoir of particles, resulting in a Poisson distribution of the variable 
z/j. This distribution is characterized by its mean value Vi, which, for reasons 
of consistency with the deterministic version, should be proportional to the 
weight a^' (see Eq. [TTj) . Normalization requires 

p. = (21) 
A* 

Equation [20] then results from assuming that all the velocity bins on the site 
are statistically independent, except for the constraints of conserved mass 
and momentum, which are taken into account by the delta functions, in 
close analogy to the statistical description of the microcanonical ensemble 
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The further development is somewhat technical but straightforward and 
shall be sketched only briefly. We use Stirling's formula and transform back 
to the rii to write the factor in front of the delta functions as exp(5'), where 
the entropy S has a Boltzmann-like form. Maximizing P is equivalent to 
maximizing S under the constraints of given values for p and j, and the 
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solution of this problem, up to second order in u, is just Eq. [TTl as is well- 
known from previous studies of the "entropic lattice Boltzmann" approach 



12l |. Fluctuations around the most probable populations are described by 
n^'^'^, which, within a saddle-point approximation, obey a Gaussian distri- 
bution, whose variance is, within a. u —>■ approximation, given by fipa'^'. 
Normalizing the fluctuations to unit variance, followed by an orthonormal 
transformation to normalized modes tti^'^'^, yields a very simple form for the 
probability distribution, 

P({m,"^n) ocexp ( -i 5^ mrM, (22) 

V fe>3 / 

where modes i = 0, . . . , 3 do not occur due to mass and momentum conser- 
vation. These modes are updated according to the rule 

* *neq ^ neq , /r)o\ 

rrif, = jkirii^ + (pkrk, (23) 

with adjustable parameters 7^, (pk, and normalized, independent Gaussian 
random numbers r^. It is then straightforward to show p, [sj that detailed 
balance holds exactly for 

^^={l-7lf\ (24) 

which turns out to be identical to the prescription of Adhikari et al. [sj. 
This shows that the stochastic analog of projecting out the kinetic modes 
is to sample them from scratch, and explains the non-trivial prefactor in 
the fluctuating stresses in a straightforward way. Furthermore 0, @], one 
may apply the Chapman-Enskog procedure to the stochastic version of the 
algorithm. This shows in a particularly concise way that the behavior in the 
hydrodynamic limit is given by Landau-Lifshitz fluctuating hydrodynamics 
[qJ, and that the details of the dynamics of the kinetic modes are indeed 
immaterial for the behavior in that limit, as already anticipated in Refs. 
0, [il. For practical simulations, however, one should prefer the more recent 
version which does satisfy detailed balance on the local scale as well. We 
believe that this is really an improvement that outweighs the computational 
costs, which are unfortunately not completely negligible. While simple LB 
algorithms have so few operations per collision step that they are typically 



limited by the bandwith of memory access in the streaming step [13|, this 



does not seem to be true here, where the generation of random numbers 
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L 


stresses-only 


full tliermalization 


10 


0.74 


0.47 


20 


0.66 


0.44 


30 


0.52 


0.39 



Table 1: Performance of the stochastic D3Q19 algorithm, using an implementation on a 
64-bit AMD Athlon 3500+ processor with 2.2 GHz CPU speed and 512 kB cache size. 
The program is part of the Mainz ESPResSo package. Simulations were run on 
simple-cubic lattices of size for 10^ lattice sweeps, and Gaussian random numbers were 
generated by the Box-Muller [l^ method. Performance data are given in MLUPS (million 
lattice-site updates per second). 

combined with the linear transformation to mode space and back contributes 
noticeably. In practice, one may say that the additional thermalization of the 
kinetic modes will slow down the algorithm by roughly 20% . . . 40% — at least 
this is what we observed for our D3Q19 implementation, see Table [H For 
large lattices the memory bottlenecks become more important than for small 
ones; for this reason, the simulations become systematically slower, while 
the performance difference between "stresses-only" vs. full thermalization 
becomes less pronounced. 

So far, only the case of an isothermal ideal gas has been thoroughly un- 
derstood. For the future, it is hoped that the present theoretical approach 
will also help develop an improved understanding of systems with non-trivial 
equations of state, and systems where thermal conduction and energy con- 
servation are taken into account. 
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